Illusion Game Reliability

Illusion Task (Session 1)

library(tidyverse)
library(ggdist)
library(ggside)
library(easystats)
library(patchwork)
library(brms)

logmod <- function(x) sign(x) * log(1 + abs(x))
sqrtmod <- function(x) sign(x) * sqrt(abs(x))
cbrtmod <- function(x) sign(x) * (abs(x)**(1 / 3))

df <- read.csv("../data/preprocessed_illusion1.csv") |>
  mutate(
    Block = as.factor(Block),
    Illusion_Side = as.factor(Illusion_Side)
  )

for(ill in c("Ebbinghaus", 
             "VerticalHorizontal", 
             "MullerLyer")) {
  print(ill)
  model <- brms::brm(
    Error ~ t2(Illusion_Difference, Illusion_Strength, k=8, bs="cr") + (1 | Participant),
    data = filter(df, Illusion_Type == ill),
    family = "bernoulli",
    algorithm="sampling"
  )

  name <- paste0("gam_", tolower(ill), "_err")
  assign(name, model)  # Rename with string
  save(list=name, file=paste0("models/", name, ".Rdata"))
  
  model <- brms::brm(
    brms::bf(
      RT ~ t2(Illusion_Difference, Illusion_Strength, k=8, bs="cr") + (1 | Participant)
      # sigma ~ t2(Illusion_Difference, Illusion_Strength, k=8, bs="cr") + poly(ISI, 2) + (1 | Participant),
      # beta ~ t2(Illusion_Difference, Illusion_Strength, k=8, bs="cr") + poly(ISI, 2) + (1 | Participant)
    ),
    data = filter(df, Illusion_Type == ill, Error == 0),
    # family = "exgaussian",
    # init=0,
    algorithm="sampling"
  )
  
  name <- paste0("gam_", tolower(ill), "_rt")
  assign(name, model)  # Rename with string
  save(list=name, file=paste0("models/", name, ".Rdata"))
}

Ebbinghaus

Error Rate

data <- filter(df, Illusion_Type == "Ebbinghaus") 

Descriptive

plot_desc_errors <- function(data) {
  dat <- data |> 
    mutate(Illusion_Difference = 
      datawizard::categorize(Illusion_Difference,
                             n_groups = 3, 
                             split = "equal_length",
                             label = "median"))
  
  col <- colorRampPalette(c("#F44336", "#FFC107", "#4CAF50"))(length(unique(dat$Illusion_Difference)))
  
  dat |> 
    ggplot(aes(x = Illusion_Strength)) +
    geom_histogram(data=filter(dat, Error == 1),
                   aes(y=..count../sum(..count..), fill = Illusion_Difference),
                   binwidth = diff(range(data$Illusion_Strength)) / 10,
                   alpha = 1/3) +
    geom_smooth(aes(y = Error, color = Illusion_Difference, fill=Illusion_Difference), 
                method = 'gam', 
                formula = y ~ s(x, bs = "cr"),
                method.args = list(family = "binomial"), 
                alpha = 1/3) +
    scale_x_continuous(expand = c(0, 0)) +
    scale_y_continuous(expand = c(0, 0), labels = scales::percent) +
    scale_color_manual(values = col) +
    scale_fill_manual(values = col) +
    coord_cartesian(ylim = c(0, 1), xlim = range(dat$Illusion_Strength)) +
    labs(x = "Illusion Strength", y = "Probability of Error") +
    theme_modern() +
    facet_grid(~Illusion_Side, labeller = "label_both")
}

plot_desc_errors(data)

Model Selection

test_models <- function(data, y = "RT") {
  # TODO: add random effect
  models <- list()
  for(diff in c("Illusion_Difference", 
             "logmod(Illusion_Difference)", 
             "sqrtmod(Illusion_Difference)", 
             "cbrtmod(Illusion_Difference)")) {
    for(ill in c("abs(Illusion_Strength)",
                 "logmod(abs(Illusion_Strength))",
                 "sqrtmod(abs(Illusion_Strength))",
                 "cbrtmod(abs(Illusion_Strength))")) {
      m <- glmmTMB::glmmTMB(
        as.formula(
          paste0(y, 
                 " ~ Illusion_Effect / (", 
                 diff,
                 " * ",
                 ill,
                 ") + (1|Participant)")),
        data = data, 
        family = ifelse(y == "RT", "gaussian", "binomial")
          )
      models[[paste(diff, "*", ill)]] <- m
    }
  }
  performance::test_performance(models) |> 
    data_select(-Model) |> 
    arrange(desc(BF))
}

insight::display(test_models(data, "Error"))

Name	BF
cbrtmod(Illusion_Difference) * abs(Illusion_Strength)	2.88
sqrtmod(Illusion_Difference) * abs(Illusion_Strength)	2.30
logmod(Illusion_Difference) * abs(Illusion_Strength)	1.30
cbrtmod(Illusion_Difference) * logmod(abs(Illusion_Strength))	0.672
sqrtmod(Illusion_Difference) * logmod(abs(Illusion_Strength))	0.536
logmod(Illusion_Difference) * logmod(abs(Illusion_Strength))	0.303
cbrtmod(Illusion_Difference) * sqrtmod(abs(Illusion_Strength))	0.270
Illusion_Difference * logmod(abs(Illusion_Strength))	0.234
sqrtmod(Illusion_Difference) * sqrtmod(abs(Illusion_Strength))	0.216
logmod(Illusion_Difference) * sqrtmod(abs(Illusion_Strength))	0.123
cbrtmod(Illusion_Difference) * cbrtmod(abs(Illusion_Strength))	0.099
Illusion_Difference * sqrtmod(abs(Illusion_Strength))	0.095
sqrtmod(Illusion_Difference) * cbrtmod(abs(Illusion_Strength))	0.079
logmod(Illusion_Difference) * cbrtmod(abs(Illusion_Strength))	0.045
Illusion_Difference * cbrtmod(abs(Illusion_Strength))	0.035
Illusion_Difference * abs(Illusion_Strength)

Each model is compared to cbrtmod(Illusion_Difference) * abs(Illusion_Strength).

Model Specification

formula <- brms::bf(
  Error ~ Illusion_Effect / (cbrtmod(Illusion_Difference) * abs(Illusion_Strength)) +
    (cbrtmod(Illusion_Difference) * abs(Illusion_Strength) | Participant),
  family = "bernoulli"
)

# brms::get_prior(formula, data = data)
# brms::validate_prior(formula)

illusion1_ebbinghaus_err <- brms::brm(formula,
  data = data,
  refresh = 100
)

save(illusion1_ebbinghaus_err, file="models/illusion1_ebbinghaus_err.Rdata")

Model Inspection

load("models/illusion1_ebbinghaus_err.Rdata")
load("models/gam_ebbinghaus_err.Rdata")

plot_model_errors <- function(data, model, gam) {
  pred <- estimate_relation(model, length = c(3, 100))
  
  dat <- data |> 
    mutate(Illusion_Difference = datawizard::categorize(
      Illusion_Difference,
      n_groups = 3, 
      split = "equal_length",
      label = "median"),
      Illusion_Strength = as.character(Illusion_Strength)) |> 
    group_by(Illusion_Strength, Illusion_Difference) |>
    summarize(Error = sum(Error),
              n = n()) |> 
    group_by(Illusion_Strength) |> 
    mutate(n = sum(n), 
           Error = Error / n,
           Illusion_Strength = as.numeric(Illusion_Strength)) |> 
    ungroup()
  
  pred |> 
    ggplot(aes(x = Illusion_Strength, y = Predicted)) +
    geom_bar(data=dat, aes(y=Error, fill = Illusion_Difference), 
             stat="identity", position=position_dodge2(reverse=TRUE), alpha = 0.2) +
    scale_fill_manual(values = c("#F44336", "#FFC107", "#4CAF50"), guide = "none") +
    ggnewscale::new_scale_fill() +
    geom_ribbon(aes(ymin = CI_low, ymax = CI_high, 
                    fill = Illusion_Difference, group = Illusion_Difference),
                alpha = 0.2) +
    geom_vline(xintercept = 0, linetype = "dashed") +
    geom_hline(yintercept = c(0.5), linetype = "dotted", alpha = 0.5) +
    geom_line(data = estimate_relation(gam, length = c(3, 100)),
              aes(color = Illusion_Difference, group = Illusion_Difference),
              linetype = "dotted") +
    geom_line(aes(color = Illusion_Difference, group = Illusion_Difference)) +
    scale_y_continuous(limits = c(0, 1), expand = c(0, 0), labels = scales::percent) +
    scale_x_continuous(expand = c(0, 0)) +
    scale_color_gradientn(colours = c("#4CAF50", "#FFC107", "#F44336"), trans = 'reverse') +
    scale_fill_gradientn(colours = c("#4CAF50", "#FFC107", "#F44336"), trans = 'reverse') +
    theme_modern(axis.title.space = 5) +
    guides(color = "none") +
    labs(
      color = "Difficulty",
      fill = "Difficulty",
      y = "Probability of Error",
      x = "Illusion Strength"
    ) +
    theme(plot.title = element_text(face = "bold", hjust = 0.5))
}

p_ebbinghaus_err <- plot_model_errors(data, illusion1_ebbinghaus_err, gam_ebbinghaus_err)
p_ebbinghaus_err

Model Performance

performance::performance(illusion1_ebbinghaus_err, metrics = c("R2", "ICC")) |> 
  display()

R2	R2 (marg.)	ICC
0.21	0.15	0.19

Response Time

data <- filter(df, Illusion_Type == "Ebbinghaus", Error == 0) 

Descriptive

plot_desc_rt <- function(data) {
  dat <- data |> 
    mutate(Illusion_Difference = 
      datawizard::categorize(Illusion_Difference,
                             n_groups = 3, 
                             split = "equal_length",
                             label = "median")) 
  
  dat |> 
    ggplot(aes(x = Illusion_Strength, y = RT)) +
    geom_jitter2(aes(color = Illusion_Difference), width = diff(range(dat$Illusion_Strength)) / 50) +
    geom_smooth(aes(y = RT, color = Illusion_Difference, fill=Illusion_Difference), 
                method = 'gam', 
                formula = y ~ s(x, bs = "cr"), 
                alpha = 1/3) +
    scale_x_continuous(expand = c(0, 0)) +
    scale_y_continuous(expand = c(0, 0)) +
    scale_color_manual(values = c("#F44336", "#FFC107", "#4CAF50")) +
    scale_fill_manual(values = c("#F44336", "#FFC107", "#4CAF50")) +
    labs(x = "Illusion Strength", y = "Response Time") +
    theme_modern() +
    ggside::geom_ysidedensity(aes(fill = Illusion_Difference), color = NA, alpha = 0.3) +
    ggside::theme_ggside_void() +
    ggside::scale_ysidex_continuous(expand = c(0, 0)) +
    ggside::ggside() +
    facet_grid(~Illusion_Side, labeller = "label_both")
}

plot_desc_rt(data)

Model Selection

insight::display(test_models(data, "RT"))

Name	BF
cbrtmod(Illusion_Difference) * abs(Illusion_Strength)	1.31
sqrtmod(Illusion_Difference) * abs(Illusion_Strength)	1.24
cbrtmod(Illusion_Difference) * logmod(abs(Illusion_Strength))	1.06
logmod(Illusion_Difference) * abs(Illusion_Strength)	1.05
sqrtmod(Illusion_Difference) * logmod(abs(Illusion_Strength))	1.01
cbrtmod(Illusion_Difference) * sqrtmod(abs(Illusion_Strength))	0.966
sqrtmod(Illusion_Difference) * sqrtmod(abs(Illusion_Strength))	0.925
logmod(Illusion_Difference) * logmod(abs(Illusion_Strength))	0.868
cbrtmod(Illusion_Difference) * cbrtmod(abs(Illusion_Strength))	0.861
Illusion_Difference * logmod(abs(Illusion_Strength))	0.835
sqrtmod(Illusion_Difference) * cbrtmod(abs(Illusion_Strength))	0.828
logmod(Illusion_Difference) * sqrtmod(abs(Illusion_Strength))	0.803
Illusion_Difference * sqrtmod(abs(Illusion_Strength))	0.776
logmod(Illusion_Difference) * cbrtmod(abs(Illusion_Strength))	0.726
Illusion_Difference * cbrtmod(abs(Illusion_Strength))	0.705
Illusion_Difference * abs(Illusion_Strength)

Each model is compared to cbrtmod(Illusion_Difference) * abs(Illusion_Strength).

Model Specification

# TODO: shift to lognormal
formula <- brms::bf(
  RT ~ Illusion_Effect / (cbrtmod(Illusion_Difference) * abs(Illusion_Strength)) +
    (Illusion_Effect / (cbrtmod(Illusion_Difference) * abs(Illusion_Strength)) | Participant)
)

# brms::get_prior(formula, data = data)
# brms::validate_prior(formula)

illusion1_ebbinghaus_rt <- brms::brm(formula,
  data = data,
  refresh = 0
)

save(illusion1_ebbinghaus_rt, file="models/illusion1_ebbinghaus_rt.Rdata")

Model Inspection

load("models/illusion1_ebbinghaus_rt.Rdata")
load("models/gam_ebbinghaus_rt.Rdata")

plot_model_rt <- function(data, model, gam) {
  pred <- estimate_relation(model, length = c(3, 100)) 
  
  dat <- data |> 
    mutate(Illusion_Difference = datawizard::categorize(
      Illusion_Difference,
      n_groups = 3, 
      split = "equal_length",
      label = "median"))
  
  col <- colorRampPalette(c("#F44336", "#FFC107", "#4CAF50"))(length(unique(dat$Illusion_Difference)))
  
  pred |> 
    ggplot(aes(x = Illusion_Strength, y = Predicted)) +
    ggdist::stat_gradientinterval(
      data=dat, 
      aes(y = RT, fill = Illusion_Difference),
      geom = "slab",
      scale = 1,
      width = diff(range(dat$Illusion_Strength)) / 20,
      fill_type = "gradient",
      position = "dodge"
    ) +
    scale_fill_manual(values = col, guide = "none") +
    ggnewscale::new_scale_fill() +
    geom_ribbon(aes(ymin = CI_low, ymax = CI_high, 
                    fill = Illusion_Difference, group = Illusion_Difference),
                alpha = 0.2) +
    geom_vline(xintercept = 0, linetype = "dashed") +
    geom_line(aes(color = Illusion_Difference, group = Illusion_Difference)) +
    geom_line(data = estimate_relation(gam, length = c(3, 100)),
              aes(color = Illusion_Difference, group = Illusion_Difference),
              linetype = "dotted") +
    scale_y_continuous(expand = c(0, 0)) +
    scale_x_continuous(expand = c(0, 0)) +
    scale_color_gradientn(colours = rev(col), trans = 'reverse') +
    scale_fill_gradientn(colours = rev(col), trans = 'reverse') +
    labs(x = "Illusion Strength", 
         y = "Response Time (s)",
         color = "Difficulty",
         fill = "Difficulty") +
    theme_modern(axis.title.space = 5) 
    # ggnewscale::new_scale_fill() +
    # scale_fill_manual(values = c("#F44336", "#FFC107", "#4CAF50"), guide = "none") +
    # ggside::geom_ysidedensity(data=dat,
    #                           aes(fill = Illusion_Difference, y = RT), color = NA, alpha = 0.3) +
    # ggside::theme_ggside_void() +
    # ggside::scale_ysidex_continuous(expand = c(0, 0)) +
    # ggside::ggside()
}

p_ebbinghaus_rt <- plot_model_rt(data, illusion1_ebbinghaus_rt, gam_ebbinghaus_rt)
p_ebbinghaus_rt

Model Performance

performance::performance(illusion1_ebbinghaus_rt, metrics = c("R2", "ICC")) |> 
  display()

R2	R2 (marg.)	ICC
0.26	0.05	0.33

Visualization

plot_all <- function(data, p_err, p_rt) {
  illname <- unique(data$Illusion_Type)
  # Get stimuli
  dat <- data |>
    filter(Error == 0) |>
    filter(
      Illusion_Type == illname,
      Answer %in% c("left", "up")
    ) |>
    select(Stimulus, Illusion_Strength, Illusion_Difference)

  dat <- rbind(
    filter(dat, Illusion_Strength == min(Illusion_Strength)) |>
      filter(Illusion_Difference %in% c(min(Illusion_Difference), max(Illusion_Difference))),
    filter(dat, Illusion_Strength == max(Illusion_Strength)) |>
      filter(Illusion_Difference %in% c(min(Illusion_Difference), max(Illusion_Difference))),
    filter(dat, Illusion_Difference == min(Illusion_Difference)) |>
      filter(Illusion_Strength %in% c(min(Illusion_Strength), max(Illusion_Strength))),
    filter(dat, Illusion_Difference == max(Illusion_Difference)) |>
      filter(Illusion_Strength %in% c(min(Illusion_Strength), max(Illusion_Strength)))
  )


  img_leftdown <- filter(dat, Illusion_Difference == max(Illusion_Difference)) |>
    filter(Illusion_Strength == min(Illusion_Strength)) |>
    pull(Stimulus) |>
    unique()
  img_rightdown <- filter(dat, Illusion_Difference == max(Illusion_Difference)) |>
    filter(Illusion_Strength == max(Illusion_Strength)) |>
    pull(Stimulus) |>
    unique()
  img_leftup <- filter(dat, Illusion_Strength == min(Illusion_Strength)) |>
    filter(Illusion_Difference == min(Illusion_Difference)) |>
    pull(Stimulus) |>
    unique()
  img_rightup <- filter(dat, Illusion_Strength == max(Illusion_Strength)) |>
    filter(Illusion_Difference == min(Illusion_Difference)) |>
    pull(Stimulus) |>
    unique()


  img_leftdown <- paste0("../session1/stimuli/", img_leftdown, ".png") |>
    png::readPNG() |>
    grid::rasterGrob(interpolate = TRUE) |>
    patchwork::wrap_elements()
  img_rightdown <- paste0("../session1/stimuli/", img_rightdown, ".png") |>
    png::readPNG() |>
    grid::rasterGrob(interpolate = TRUE) |>
    patchwork::wrap_elements()
  img_leftup <- paste0("../session1/stimuli/", img_leftup, ".png") |>
    png::readPNG() |>
    grid::rasterGrob(interpolate = TRUE) |>
    patchwork::wrap_elements()
  img_rightup <- paste0("../session1/stimuli/", img_rightup, ".png") |>
    png::readPNG() |>
    grid::rasterGrob(interpolate = TRUE) |>
    patchwork::wrap_elements()

  p <- ((p_err + theme(axis.title.x = element_blank(),
                       plot.title = element_blank())) /
          (p_rt + theme(plot.title = element_blank(), legend.position = "none"))) +
    patchwork::plot_layout(guides = "collect")

  wrap_elements(((img_leftup | patchwork::plot_spacer() | img_rightup) / p / (img_leftdown | patchwork::plot_spacer() | img_rightdown) +
    patchwork::plot_layout(heights = c(0.5, 1.5, 0.5)) +
    patchwork::plot_annotation(title = case_when(
      illname == "MullerLyer" ~ "Müller-Lyer",
      illname == "VerticalHorizontal" ~ "Vertical-Horizontal",
      TRUE ~ "Ebbinghaus"
      ),
      theme = theme(plot.title = element_text(size=rel(1.75), face="bold", hjust=0.5)))))
}

p_ebbinghaus <- plot_all(data, p_ebbinghaus_err, p_ebbinghaus_rt)
p_ebbinghaus

Müller-Lyer

Error Rate

data <- filter(df, Illusion_Type == "MullerLyer") 

Descriptive

plot_desc_errors(data)

Model Selection

insight::display(test_models(data, "Error"))

Name	BF
Illusion_Difference * logmod(abs(Illusion_Strength))	47.80
Illusion_Difference * cbrtmod(abs(Illusion_Strength))	33.77
logmod(Illusion_Difference) * logmod(abs(Illusion_Strength))	29.55
logmod(Illusion_Difference) * cbrtmod(abs(Illusion_Strength))	21.99
Illusion_Difference * sqrtmod(abs(Illusion_Strength))	19.38
logmod(Illusion_Difference) * sqrtmod(abs(Illusion_Strength))	13.06
sqrtmod(Illusion_Difference) * logmod(abs(Illusion_Strength))	1.71
sqrtmod(Illusion_Difference) * cbrtmod(abs(Illusion_Strength))	1.49
sqrtmod(Illusion_Difference) * sqrtmod(abs(Illusion_Strength))	0.979
logmod(Illusion_Difference) * abs(Illusion_Strength)	0.728
cbrtmod(Illusion_Difference) * logmod(abs(Illusion_Strength))	0.261
cbrtmod(Illusion_Difference) * cbrtmod(abs(Illusion_Strength))	0.241
cbrtmod(Illusion_Difference) * sqrtmod(abs(Illusion_Strength))	0.166
sqrtmod(Illusion_Difference) * abs(Illusion_Strength)	0.069
cbrtmod(Illusion_Difference) * abs(Illusion_Strength)	0.013
Illusion_Difference * abs(Illusion_Strength)

Each model is compared to Illusion_Difference * logmod(abs(Illusion_Strength)).

Model Specification

formula <- brms::bf(
  Error ~ Illusion_Effect / (Illusion_Difference * logmod(abs(Illusion_Strength))) +
     (Illusion_Difference * logmod(abs(Illusion_Strength)) | Participant),
  family = "bernoulli"
)

# brms::get_prior(formula, data = data)
# brms::validate_prior(formula)

illusion1_mullerlyer_err <- brms::brm(formula,
  data = data,
  refresh = 100
)

save(illusion1_mullerlyer_err, file="models/illusion1_mullerlyer_err.Rdata")

Model Inspection

load("models/illusion1_mullerlyer_err.Rdata")
load("models/gam_mullerlyer_err.Rdata")

p_mullerlyer_err <- plot_model_errors(data, illusion1_mullerlyer_err, gam_mullerlyer_err)
p_mullerlyer_err

Model Performance

performance::performance(illusion1_mullerlyer_err, metrics = c("R2", "ICC")) |> 
  display()

R2	R2 (marg.)	ICC
0.69	0.64	0.29

Response Time

data <- filter(df, Illusion_Type == "MullerLyer", Error == 0) 

Descriptive

plot_desc_rt(data)

Model Selection

insight::display(test_models(data, "RT"))

Name	BF
cbrtmod(Illusion_Difference) * abs(Illusion_Strength)	61.55
cbrtmod(Illusion_Difference) * sqrtmod(abs(Illusion_Strength))	40.70
sqrtmod(Illusion_Difference) * abs(Illusion_Strength)	28.29
sqrtmod(Illusion_Difference) * sqrtmod(abs(Illusion_Strength))	17.74
cbrtmod(Illusion_Difference) * cbrtmod(abs(Illusion_Strength))	15.88
sqrtmod(Illusion_Difference) * cbrtmod(abs(Illusion_Strength))	6.83
cbrtmod(Illusion_Difference) * logmod(abs(Illusion_Strength))	3.31
logmod(Illusion_Difference) * abs(Illusion_Strength)	2.87
logmod(Illusion_Difference) * sqrtmod(abs(Illusion_Strength))	1.63
sqrtmod(Illusion_Difference) * logmod(abs(Illusion_Strength))	1.41
logmod(Illusion_Difference) * cbrtmod(abs(Illusion_Strength))	0.620
Illusion_Difference * sqrtmod(abs(Illusion_Strength))	0.551
Illusion_Difference * cbrtmod(abs(Illusion_Strength))	0.209
logmod(Illusion_Difference) * logmod(abs(Illusion_Strength))	0.126
Illusion_Difference * logmod(abs(Illusion_Strength))	0.042
Illusion_Difference * abs(Illusion_Strength)

Each model is compared to cbrtmod(Illusion_Difference) * abs(Illusion_Strength).

Model Specification

# TODO: shift to lognormal
formula <- brms::bf(
  RT ~ Illusion_Effect / (cbrtmod(Illusion_Difference) * abs(Illusion_Strength)) +
    (Illusion_Effect / (cbrtmod(Illusion_Difference) * abs(Illusion_Strength)) | Participant)
)

# brms::get_prior(formula, data = data)
# brms::validate_prior(formula)

illusion1_mullerlyer_rt <- brms::brm(formula,
  data = data,
  refresh = 0
)

save(illusion1_mullerlyer_rt, file="models/illusion1_mullerlyer_rt.Rdata")

Model Inspection

load("models/illusion1_mullerlyer_rt.Rdata")
load("models/gam_mullerlyer_rt.Rdata")

p_mullerlyer_rt <- plot_model_rt(data, illusion1_mullerlyer_rt, gam_mullerlyer_rt)
p_mullerlyer_rt

Model Performance

performance::performance(illusion1_mullerlyer_rt, metrics = c("R2", "ICC")) |> 
  display()

R2	R2 (marg.)	ICC
0.45	0.29	0.38

Visualization

p_mullerlyer <- plot_all(data, p_mullerlyer_err, p_mullerlyer_rt)
p_mullerlyer

Vertical-Horizontal

Error Rate

data <- filter(df, Illusion_Type == "VerticalHorizontal") 

Descriptive

plot_desc_errors(data)

Model Selection

insight::display(test_models(data, "Error"))

Name	BF
sqrtmod(Illusion_Difference) * sqrtmod(abs(Illusion_Strength))	1.91
cbrtmod(Illusion_Difference) * sqrtmod(abs(Illusion_Strength))	1.90
sqrtmod(Illusion_Difference) * abs(Illusion_Strength)	1.63
cbrtmod(Illusion_Difference) * abs(Illusion_Strength)	1.63
sqrtmod(Illusion_Difference) * cbrtmod(abs(Illusion_Strength))	1.61
cbrtmod(Illusion_Difference) * cbrtmod(abs(Illusion_Strength))	1.60
logmod(Illusion_Difference) * sqrtmod(abs(Illusion_Strength))	1.43
logmod(Illusion_Difference) * cbrtmod(abs(Illusion_Strength))	1.21
logmod(Illusion_Difference) * abs(Illusion_Strength)	1.21
Illusion_Difference * sqrtmod(abs(Illusion_Strength))	1.19
Illusion_Difference * cbrtmod(abs(Illusion_Strength))	1.01
sqrtmod(Illusion_Difference) * logmod(abs(Illusion_Strength))	0.911
cbrtmod(Illusion_Difference) * logmod(abs(Illusion_Strength))	0.902
logmod(Illusion_Difference) * logmod(abs(Illusion_Strength))	0.691
Illusion_Difference * logmod(abs(Illusion_Strength))	0.576
Illusion_Difference * abs(Illusion_Strength)

Each model is compared to sqrtmod(Illusion_Difference) * sqrtmod(abs(Illusion_Strength)).

Model Specification

formula <- brms::bf(
  Error ~ Illusion_Effect / (sqrtmod(Illusion_Difference) * sqrtmod(abs(Illusion_Strength))) +
     (sqrtmod(Illusion_Difference) * sqrtmod(abs(Illusion_Strength)) | Participant),
  family = "bernoulli"
)

# brms::get_prior(formula, data = data)
# brms::validate_prior(formula)

illusion1_verticalhorizontal_err <- brms::brm(formula,
  data = data,
  refresh = 100
)

save(illusion1_verticalhorizontal_err, file="models/illusion1_verticalhorizontal_err.Rdata")

Model Inspection

load("models/illusion1_verticalhorizontal_err.Rdata")
load("models/gam_verticalhorizontal_err.Rdata")

p_verticalhorizontal_err <- plot_model_errors(data, 
                                              illusion1_verticalhorizontal_err, 
                                              gam_verticalhorizontal_err)
p_verticalhorizontal_err

Model Performance

performance::performance(illusion1_verticalhorizontal_err, metrics = c("R2", "ICC")) |> 
  display()

R2	R2 (marg.)	ICC
0.30	0.24	0.17

Response Time

data <- filter(df, Illusion_Type == "VerticalHorizontal", Error == 0) 

Descriptive

plot_desc_rt(data)

Model Selection

insight::display(test_models(data, "RT"))

Name	BF
cbrtmod(Illusion_Difference) * logmod(abs(Illusion_Strength))	10.03
cbrtmod(Illusion_Difference) * cbrtmod(abs(Illusion_Strength))	9.56
sqrtmod(Illusion_Difference) * logmod(abs(Illusion_Strength))	8.80
sqrtmod(Illusion_Difference) * cbrtmod(abs(Illusion_Strength))	8.40
cbrtmod(Illusion_Difference) * sqrtmod(abs(Illusion_Strength))	7.99
sqrtmod(Illusion_Difference) * sqrtmod(abs(Illusion_Strength))	7.03
logmod(Illusion_Difference) * logmod(abs(Illusion_Strength))	5.05
logmod(Illusion_Difference) * cbrtmod(abs(Illusion_Strength))	4.82
Illusion_Difference * logmod(abs(Illusion_Strength))	4.08
logmod(Illusion_Difference) * sqrtmod(abs(Illusion_Strength))	4.05
Illusion_Difference * cbrtmod(abs(Illusion_Strength))	3.90
Illusion_Difference * sqrtmod(abs(Illusion_Strength))	3.27
cbrtmod(Illusion_Difference) * abs(Illusion_Strength)	2.41
sqrtmod(Illusion_Difference) * abs(Illusion_Strength)	2.13
logmod(Illusion_Difference) * abs(Illusion_Strength)	1.23
Illusion_Difference * abs(Illusion_Strength)

Each model is compared to cbrtmod(Illusion_Difference) * logmod(abs(Illusion_Strength)).

Model Specification

# TODO: shift to lognormal
formula <- brms::bf(
  RT ~ Illusion_Effect / (cbrtmod(Illusion_Difference) * cbrtmod(abs(Illusion_Strength))) +
    (Illusion_Effect / (cbrtmod(Illusion_Difference) * cbrtmod(abs(Illusion_Strength))) | Participant)
)

# brms::get_prior(formula, data = data)
# brms::validate_prior(formula)

illusion1_verticalhorizontal_rt <- brms::brm(formula,
  data = data,
  refresh = 0
)

save(illusion1_verticalhorizontal_rt, file="models/illusion1_verticalhorizontal_rt.Rdata")

Model Inspection

load("models/illusion1_verticalhorizontal_rt.Rdata")
load("models/gam_verticalhorizontal_rt.Rdata")

p_verticalhorizontal_rt <- plot_model_rt(data, 
                                         illusion1_verticalhorizontal_rt, 
                                         gam_verticalhorizontal_rt)
p_verticalhorizontal_rt

Model Performance

performance::performance(illusion1_verticalhorizontal_rt, metrics = c("R2", "ICC")) |> 
  display()

R2	R2 (marg.)	ICC
0.46	0.17	0.61

Visualization

p_verticalhorizontal <- plot_all(data, p_verticalhorizontal_err, p_verticalhorizontal_rt)
p_verticalhorizontal

Individual Scores

p <- p_ebbinghaus | p_mullerlyer | p_verticalhorizontal 
p

# ggsave("figures/figure1.png", p, dpi=300, width=16, height=10)

get_scores <- function(model, illusion="Ebbinghaus") {
  family <- insight::find_response(model)
  scores <- modelbased::estimate_grouplevel(model) |> 
    data_filter(str_detect(Level, "Participant")) |> 
    data_filter(!str_detect(Group, "Illusion_EffectCongruent|Intercept")) |> 
    mutate(Group = str_remove_all(Group, ": Participant"),
           Level = str_remove_all(Level, "Participant."),
           Group = str_remove_all(Group, "Illusion_EffectIncongruent:"),
           Group = str_remove_all(Group, "cbrtmod"),
           Group = str_remove_all(Group, "sqrtmod"),
           Group = str_remove_all(Group, "logmod"),
           Group = str_remove_all(Group, "abs"),
           Group = str_replace(Group, "Illusion_Difference:Illusion_Strength", "Interaction"),
           Group = str_replace(Group, "Illusion_Difference", "Diff"),
           Group = str_replace(Group, "Illusion_Strength", "Strength")) |>
    data_filter(!str_detect(Group, "Illusion_EffectIncongruent"))
  
  p <- scores |> 
    ggplot(aes(x = Median, y = Level)) +
    geom_pointrange(aes(xmin = CI_low, xmax = CI_high, color = Group)) +
    scale_color_flat_d() +
    scale_fill_flat_d() +
    labs(y = "Participants") +
    theme_modern() +
    theme(strip.placement = "oustide",
          axis.title.x = element_blank(),
          axis.text.y = element_blank()) +
    ggside::geom_xsidedensity(aes(fill=Group, y = after_stat(scaled)), color = NA, alpha = 0.3) +
    ggside::theme_ggside_void() +
    ggside::scale_xsidey_continuous(expand = c(0, 0)) +
    ggside::ggside() +
    facet_grid(~Group, switch = "both", scales = "free")  +
    ggtitle(paste(illusion, "-", family))
  
  # Reshape
  scores <- scores |> 
    select(Group, Participant = Level, Median) |> 
    pivot_wider(names_from = "Group", values_from = "Median") |> 
    data_addprefix(paste0(illusion, "_"), select=-1) |> 
    data_addsuffix(paste0("_", family), select=-1)
  list(scores = scores, p = p)
}

ebbinghaus_err <- get_scores(illusion1_ebbinghaus_err, illusion="Ebbinghaus")
ebbinghaus_rt <- get_scores(illusion1_ebbinghaus_rt, illusion="Ebbinghaus")
mullerlyer_err <- get_scores(illusion1_mullerlyer_err, illusion="MullerLyer")
mullerlyer_rt <- get_scores(illusion1_mullerlyer_rt, illusion="MullerLyer")
verticalhorizontal_err <- get_scores(illusion1_verticalhorizontal_err, illusion="VerticalHorizontal")
verticalhorizontal_rt <- get_scores(illusion1_verticalhorizontal_rt, illusion="VerticalHorizontal")

p <- (ebbinghaus_err$p + ebbinghaus_rt$p) /
  (mullerlyer_err$p + mullerlyer_rt$p) /
  (verticalhorizontal_err$p + verticalhorizontal_rt$p) +
  plot_layout(guides = "collect") +
  plot_annotation(title = "Inter- and Intra- Variability of Illusion Scores", theme = theme(plot.title = element_text(face = "bold", hjust = 0.5))) 
p

scores <- ebbinghaus_err$scores |> 
  merge(ebbinghaus_rt$scores, by="Participant") |> 
  merge(mullerlyer_err$scores, by="Participant") |> 
  merge(mullerlyer_rt$scores, by="Participant") |> 
  merge(verticalhorizontal_err$scores, by="Participant") |> 
  merge(verticalhorizontal_rt$scores, by="Participant")

write.csv(scores, "../data/scores_illusion1.csv", row.names = FALSE)